Possible depth-resolved reconstruction of shear moduli in the cornea following collagen crosslinking (CXL) with optical coherence tomography and elastography

Collagen crosslinking of the cornea (CXL) is commonly employed to prevent or treat keratoconus. Although the change of corneal stiffness induced by CXL surgery can be monitored with non-contact dynamic Optical Coherence Elastography (OCE) by tracking mechanical wave propagation, the depth dependence of this change is still unclear if the cornea is not crosslinked through the whole depth. Here we propose to combine phase-decorrelation measurement applied to OCT structural images and acoustic micro-tapping (A$\mu$T) OCE to explore possible depth reconstruction of stiffness within crosslinked corneas in an ex vivo human cornea sample. The analysis of experimental OCT images is used to define the penetration depth of CXL into the cornea, which varies from $\sim$100$\mu m$ in the periphery to $\sim$150$\mu m$ in the central area and exhibits a sharp transition between areas. This information was used in a two-layer analytical model to quantify the stiffness of the treated layer. We also discuss how the elastic moduli of partially CXL-treated cornea layers reconstructed from OCE measurements reflect the effective mechanical stiffness of the entire cornea to properly quantify surgical outcome.


I. INTRODUCTION
T the interface between air and the inner eye, the cornea provides protection and is the primary optical element focusing light onto the retina. It contains multiple layers, including epithelium and stroma. The first acts as a barrier against the external environment and the latter maintains its stiffness, transparency and focusing power [1], [2]. The microstructure of the stroma is composed of collagen fibrils, arranged in lamellae, lying within a protein rich, hydrated proteoglycan mesh [3], [4]. Corneal diseases (such as keratoconus (KC)) and surgical complications from refractive surgeries (such as LASIK) may deform the cornea (ectasia) and reduce vision. The prevalence of KC in the general population is estimated to be 1.38 per 1000 [5], and nearly 1 million refractive surgeries are performed each year in the USA. Despite their overall success, however, suboptimal visual outcomes and post-refractive corneal decompensation cannot always be predicted for an individual patient. Corneal collagen crosslinking (CXL) is a minimally invasive procedure that can potentially slow the progression of corneal ectasia [6]- [9]. UV light modifies the microstructure of the cornea soaked in riboflavin and forms additional chemical bonds between collagen fibers in the stroma [10]. Posttreatment corneas become stiffer and more resistant to enzymatic digestion [11]- [13]. Although corneal topography (curvature map) and thickness map can be obtained preoperatively, and needed refractive corrections can be estimated, there is an unmet need to predict corneal decompensation from interventions such as LASIK and CXL therapies. Unfortunately, surgical planning cannot be customized and outcomes (e.g. post-surgery corneal ectasia risks) cannot be accurately predicted without quantitatively mapping corneal elasticity. Thus, methods that can quantitatively reconstruct corneal elastic moduli are in high demand. Ocular response analyzer (ORA -Reichert Technologies) and Dynamic Scheimpflug Analyzer (DSA -Corvis ST -Oculus Opitkgerate GmbH) are the state-of-the-art in clinical measurements of corneal mechanics. They estimate stiffness as the pressure at inward applanation divided by corneal displacement [14]- [16]. Over 100 papers are published annually on the topic, and the Journal of Cataract and Refractive Surgery devoted an entire issue to the topic [15]. However, measurements induce large corneal deformations that are often clinically unacceptable, require a non-trivial IOP correction in simulations [17] and assume a simple isotropic mechanical model leading to high variability with experimental conditions. Results obtained with the Corvis ST on KC may be contradictory, and some even show no significant change in corneal stiffness pre-and post-CXL surgery [18], [19]. In addition, the result is averaged over the entire cornea with no spatial resolution, and the reconstruction is questionable if corneal thickness varies. There is no consensus in the literature on the elastic model for cornea and the stiffness range even for healthy subjects. The most common model assumes an incompressible, isotropic, and A linear elastic solid, where a single parameter, the Young's modulus (or equivalently the shear modulus ), defines elasticity. Destructive mechanical tests can determine E ex vivo, with reported values for human cornea (at low-strain) of 800 kPa to 4.7 MPa for tensile loading [12], [20]- [23], and 100 kPa to 3 MPa for inflation loading [24], [25]. Note that the destructive nature of mechanical tests precludes their clinical translation. Dynamic elastography is a promising tool to probe soft tissue biomechanics. A shear wave can be launched using direct contact excitation [26]- [29] or radiation force-based techniques [30], [31]. By tracking shear wave propagation, either using MRI [32], [33] ultrasound [26]- [28], [34] or dynamic phasesensitive OCT techniques [35]- [37], one can infer, with an appropriate mechanical description, the linear [38]- [41], or non-linear [42], [43], stiffness moduli of the tissue. Optical coherence elastography (OCE) is particularly suited to probe corneal biomechanics non-invasively in a clinical environment [30], [36], [44]- [46], as it can be combined with non-contact excitation techniques (for example, using an air-puff or acoustic micro-tapping (AμT [30])). Because cornea is thin and bounded, above by air and below by aqueous humor, wave propagation within it is guided, leading to strong geometric dispersion [39], [47]. As such, the common approach associating the Rayleigh surface wave group velocity to stiffness [30], [36], [44]- [46], [48]- [51] is not appropriate. It in fact results in a 2 order of magnitude mismatch in Young's modulus compared to tensile and inflation tests [47]. Accounting for phase-velocity dispersion is mandatory to recover corneal stiffness. Using Lamb-wave dispersion, OCE studies reported corneal shear moduli in the range of 1.8 -52.3 kPa [49], [52], [53], in close agreement with values obtained from torsional and rheometry testing of ex vivo cornea (2.5 -47.3 kPa) [54]- [56]. However, all shear-based methods produce moduli differing by 1-2 orders of magnitude from those reported by tensile and inflation tests. Recently, we hypothesized that corneal anisotropy is the primary cause of these discrepancies [39]. Corneal microstructure supports this hypothesis. The stroma contains collagen lamellae running in-plane across its width. Lamellae make up approximately 90% of tissue thickness and account for most of the cornea's mechanical structure. They are stacked vertically in approximately 200-500 separate planes [57], [58], suggesting an anisotropic mechanical behavior with very different responses to in-plane versus out-of-plane loads. We introduced a model of a nearly-incompressible transverse isotropic (NITI) medium [39], in which corneal stiffness is defined by two (in-, , and out-of-plane, ) shear moduli, decoupling tensile/inflation properties from shear responses. Based on this model, we developed an algorithm utilizing guided mechanical waves in a bounded NITI medium to reconstruct both moduli from AµT-OCE. The model was confirmed ex vivo in rabbit [59], porcine [60] and human [38], [39] models and in vivo with rabbit models [59]. In [38] we showed that both in-and out-of-plane post-CXL corneal shear moduli increased compared to their pre-surgery values, with an averaged two-fold increase in Young's modulus and an almost four-fold increase for the out-of-plane shear modulus . This confirmed that CXL better crosslinks corneal lamellae. However, deformations under physiological conditions are defined by the Young's modulus, which is less affected and has implications for potential refractive changes. Note, however, that the effects of crosslinking are always nonhomogeneous in depth, as riboflavin penetration is more pronounced in the anterior region of the cornea where the solution is applied during soaking [61], and because UVA irradiation is more efficient on the anterior part of the cornea [62]. Together, this generally leads to a clear demarcation line, indicating the efficiency of the treatment [63]. Note that riboflavin penetration is purposely limited to minimize damage to the endothelium and ensure the biological integrity of the cornea [6]. As noted above, mechanical waves generated in the cornea are guided. Such guided waves occupy the entire depth of the cornea and, therefore, carry average information over depth. Thus, reconstructing the depth dependence of corneal moduli is practically impossible without a good estimate of the penetration of CXL into the cornea. Blackburn et al. [64] have recently introduced a novel metric to track CXL penetration within the cornea using time-resolved OCT. They demonstrated that the phase decorrelation decay rate of the complex OCT signal is reduced in CXL areas and can be used to distinguish treated and untreated areas postoperatively. In this paper, we combine the methods described in [64] with AµT-OCE measurements to explore possible reconstruction of both in-and out-of-plane corneal elastic moduli over depth. For this purpose, we developed an analytical model of guided wave propagation accounting for multiple layers in the cornea, each with distinct stiffness moduli and thickness. In addition, elastic moduli prior to CXL can be directly estimated from AµT-OCE so that procedure induced changes can be quantified. We also discuss how the elastic moduli of a CXL-treated layer reconstructed from OCE measurements influences the effective mechanical stiffness of the entire cornea to predict its deformational properties and properly quantify surgical outcome.

A. Cornea preparation
One corneal-scleral ring, stored in Optisol (Chiron Ophtalmics) was obtained from CorneaGen. This sample from a 26 yearsold donor was stored for less than 30 days. The corneal-scleral button is mounted on an artificial anterior chamber (Barron, CorzaMedical; see Fig. 1), connected through the inlet port to an elevated bath filled with balanced saline solution (BSS) to apply a controlled pressure on the anterior segment of the cornea. The outlet port remained closed to allow the pressure to settle at 15 mmHg within the chamber, corresponding to human physiological conditions [65]. Crosslinking was performed following the Dresden protocol [6]. First, the epithelial membrane was removed from the sample. Then, the cornea was soaked in riboflavin for 30 min by applying a 50µL drop of 0.1% riboflavin in 20% dextran solution every two minutes. The cornea was then exposed to 3mW/cm 2 of 370 nm ultraviolet (UV) light for 30 minutes, while a drop was re-applied every 5 minutes.

B. AµT-OCE imaging system
A spectral domain OCT system with a 46.5 kHz frame rate was used to track wave propagation and structural changes within the cornea. As described in previous studies [30], [36], [38], [39], a cylindrically focused air-coupled transducer, operating at a 1 MHz frequency, generated a spatio-temporal sharp displacement at the surface of the cornea using reflection-based acoustic radiation force. Because of the transducer's cylindrical geometry, the push was line-shaped and generated quasi-planar guided waves within bounded tissue. The OCT system operated in M-B mode. A single push was triggered by the system while 512 consecutive A-scans were taken at a fixed location (Mscan). The M-scan sequence and push excitation were repeated for 256 locations, creating a three-dimensional volume with 256 -samples, 1024 -samples and 512 -samples (see Fig. 2a)), with an effective imaging range of . The vertical particle velocity was obtained from the optical phase difference between two consecutive A-lines at each location [66]. The spatio-temporal ( -) surface signature of the guided wave was computed from an exponentially weightedaverage of the vertical particle velocity over the first 180 µm of the anterior part of the cornea. As shown in Fig. 2b), the guided wave only propagated during the first 4 ms of the scans, which was used to determine the stiffness of the material by fitting the computed dispersion curve in the frequency-wavenumber domain ( -) obtained from the 2D Fourier spectrum (Fig. 2c). This procedure is detailed in Section IID. On the other hand, data from the last 7 ms were used to study structural changes with phase decorrelation (see Section IIF).

C. Multi-layer NITI model
Like most biological tissue, the cornea is nearly incompressible. In addition, its microstructure implies its transverse isotropy [60] and, therefore, its mechanical behavior under small deformation can be described with the NITI model [39]. In Voigt notation, Hook's law of stress and strain for a NITI material takes the form: , (1) where denotes engineering stress, denotes engineering strain, denotes shear stress, denotes shear strains, the subscripts , and refer to the Cartesian axes and , , and are four independent elastic constants. In previous studies [40], [60], we have demonstrated that , which accounts for tissue tensile anisotropy, cannot be determined from guided wave propagation but that the in-plane Young's modulus can be approximated as assuming tensile isotropy ( ). Thus, among the four elastic constants, only and , respectively the out-of-plane and in-plane shear moduli, are  needed to predict corneal deformation under mechanical loading. The effects of CXL on the cornea depend on depth. Several recent studies showed that postoperative CXL corneas may experience non-uniform cross-linkage with depth. The transition between crosslinked (anterior) and non-crosslinked (posterior) parts tends to be sharp rather than smooth [61][62][63]. This effect is also observed in our experiments (see Fig. 3d, e). Thus, a two-layer model, although not exactly accurate, is considered a working model to quantify postoperative corneas. CXL was also shown to change collagen fiber diameter and interfibrillar spacing [61], but nothing suggests a modification of its macroscopic anisotropic organization. Based on this observation, we developed a multi-layer model to predict wave propagation within CXL corneas (Supplementary Material 1) that accounts for any arbitrary number of layers, each with a stress-strain relationship given by Eq. (1) and linked by a perfect solid-solid boundary condition (continuity of normal components of stress and displacement across every interface). Accounting for the appropriate external boundary conditions (liquid below and air above the cornea) and the finite thickness of the medium, the dispersion relation for guided waves can be determined directly from stiffness moduli and and the thickness of each layer. A detailed description of the multi-layer model is provided in Supplementary Material 1. Although only 2 layers were considered in this study, the multilayer model will allow further refinements if, for example, the transition between crosslinked and non-crosslinked areas needs a more accurate description. Note that in an untreated cornea, only the first anti-symmetric mode, referred to as , propagates in the range of frequencies that can be recorded in elastography (typically kHz). Because a CXL cornea is approximated as two horizontally assembled layers, each having a vertically aligned symmetry axis, this symmetry holds for the global material. Based on previous experiments on CXL corneas [38], [59], and symmetry conservation, the -mode is expected to contain the primary energy of guided waves in CXL-treated tissue.

D. Fitting pre-and post-CXL
The experimental -spectrum (see Fig. 2c) was obtained by computing the 2D FFT of the -plot. The shear moduli and in pre-CXL cornea were obtained from fitting the measured -spectrum with the analytical dispersion relationship of the mode. Prior to treatment, the cornea was assumed homogeneous, which in our model corresponds to a single layer bounded above by air and below by water. Using the model described above, the dispersion curve of a single layer material was used to determine the untreated properties of the cornea. In CXL cornea, the thickness of both layers can be measured (see Section IIE) and the posterior layer can be assumed to still possess the original (i.e. untreated) elastic properties. Thus, the 2-layer model with known elastic moduli of the bottom (untreated) layer can be used to determine the stiffness of the top layer. Since the -mode carries depth averaged properties of the material, we also used a single-layer model to estimate the 'effective' stiffness of the treated cornea by fitting itsspectrum with the analytical dispersion relationship of the mode computed for a single layer. To ensure reliable fitting for all cases, we computed a goodness of fit metric ∑ ∑ , where corresponds to the energy of the 2D spectrum covered by the best analytical solution (one or two layers) at a given frequency and is the unconstrained maximal energy of the spectrum at frequency . Based on recent results (see Supplemental Material in [38]), reliable fitting in human ex vivo corneas is associated with values of . An example of a 2Dspectrum and the fitted mode obtained with this procedure for the untreated case is shown in Fig. 2c. More details about the fitting procedure, including how to determine uncertainty intervals for and for both treated and untreated cases, are given in Supplementary Material 2.

E. FEM simulations
We designed finite element (FEM) simulations in OnScale to determine the efficiency of our multi-layer NITI model for reconstruction of stiffness along depth. The geometry is shown in Fig. 4a. Corneal boundary conditions were replicated so that the material is bounded above by air and below by water. The speed of sound in all layers (material and water) was fixed to avoid reflection of compressional waves at air-tissue boundaries. It also improved the absorption of waves at the absorbing boundaries and, thus, avoided divergence of the simulations. We used transient excitation mimicking AµT experiments to generate broadband elastic waves within the material. More details about the simulations can be found in [42,50]. Based on phase-decorrelation measurements (see Section IIIF), we assumed that after CXL two layers with distinct thicknesses were formed within the cornea, the top layer being stiffer than the bottom one. The stiffness values assessed from experiments were also used in the simulations. We focused on the top surface signature of the simulated wave (see Fig. 4b) and its associated spectrum (see Fig. 4c) to study the efficiency of our fitting routine to quantify corneal elastic moduli and their variation with depth.

F. Phase Decorrelation OCT (PhD-OCT)
Blackburn et al. [64] have recently introduced a novel metric to track CXL penetration within the cornea using time-resolved OCT. It was shown that the phase decorrelation decay rate of the complex OCT signal is reduced in CXL areas and can be used to distinguish treated and untreated areas after the procedure. In our study, the autocorrelation function of the signal was computed over 15 consecutive samples at 46,500 Hz for six consecutive pixels within a given A-line: which is expected to follow an exponential decay [67]: . where is the decorrelation coefficient that is inversely proportional to the Brownian diffusion coefficient [67], meaning that the more coherent the material, the smaller the decorrelation coefficient. The procedure was performed starting at , , , … A-lines, where is the first time-sample used for phase-decorrelation ( 4 ms). The decorrelation coefficient was then computed using the averaged over the number of starting points by fitting with a first order polynomial (see Fig. 2e): , where denotes the average over the number of starting points. In crosslinked regions of the cornea (anterior), tissue stiffens, which implies that should be smaller than in the untreated region (posterior). For post-processing, we rejected all fits with b < 0.95, corresponding in general to peripheral regions where the signal to noise ratio (SNR) was too low.

A. Thickness of CXL layer
The spatial distribution of the OCT intensity signal (Figs. 3a, b, e) and phase decorrelation images (Figs. 3c, d, e) both show a clear layering effect in the treated cornea. The effect of CXL is not homogeneous across the cornea, with a more pronounced effect at the center ( ) than at the periphery ( ). For the present case, we estimated that about 30% of the cornea was treated efficiently. As shown in Table 1, the treated cornea is thinner than that prior to CXL (its thickness reduced from to ), as generally observed in the literature [68], [69].

B. Stiffness of CXL-treated corneal layer
AµT-OCE scans, taking approximatively 3 s to acquire and save data, were acquired before and after CXL. The signature of the vertical particle velocity (see Fig. 2ba) in the untreated cornea was used to compute the -spectrum (see Fig. 2cb), which was fitted using the procedure detailed in Section IIC and Supplementary Material 2 assuming the cornea as a single homogeneous layer. The fitting routine was also detailed in our recent work [38], [59]. Results for the reconstructed in-( ) and out-of-plane ( ) corneal shear moduli are shown in Table 1. It is worth mentioning that error bars are asymmetric, and the uncertainties in are always broader than those for . More details about the procedure used to compute the uncertainty intervals are given in Supplementary Material 2 (see also [38], [59]). To reconstruct the depth-dependent stiffness moduli of the cornea after it was treated with CXL, it is assumed that: i) the thickness of both layers can be measured using dynamic OCT from phase-decorrelation or intensity variation methods (see Fig. 3) which were estimated as about and with both methods; ii) the stiffness of the posterior layer is unchanged after CXL; iii) the effect of CXL is homogeneous in the anterior layer. Fixing known parameters (untreated cornea thicknesses and posterior stiffness moduli), stiffness moduli of the anterior cornea layer can be determined by fitting the wave dispersion curve in the -domain with the 2-layer model (Fig. 5b, Supplementary Material 2). The goodness of fit surface plots are shown in Fig. 4c, d. We found an increase in both the stiffness moduli and compared to that for the posterior region kPa and MPa. The goodness of fit for the 2-layer model remained high, 0.95, which is still within the range of reliable fitting. These results are in good agreement with our recent study of corneas cross-lined throughout depth [38].

C. Mixing rules of the mechanical moduli for layered materials
A theory for effective moduli of multi-layer materials was developed in the early 1970's for composite materials. It is now accepted in material science and broadly used in the development of composite structures [70]- [72]. For partially cross-linked cornea, it would be particularly interesting to have a rule to compute 'effective' moduli to estimate its ultimate stiffness, assess treatment efficiency, and predict corneal behavior and deformation post-operatively. The derivation of effective material moduli is based on 'mixing rules' of mechanical moduli across depth using the following assumptions: i) out-of-plane stresses and in-plane strains are uniform across thickness; ii) in-plane stresses and out-of-plane strains are averaged across thickness based on layer volume fractions. Note that the solution is valid only for very low frequencies. Sun et al. [71] have demonstrated that the effective low-frequency out-of-plane modulus can be computed using the inverse rule of mixture for out-of-plane material constants of individual material layers.: where is the total material thickness, is the thickness of the th layer and is the out-of-plane modulus of the th layer. On the other hand, the effective low-frequency in-plane modulus can be obtained from the rule of 'mixtures': (5) where is the in-plane modulus of the th layer. Based on the mixture rules described above, the effective lowfrequency moduli of a partially treated cornea can be computed, and in our case are equal to and .

D. Effective 'guided wave' corneal moduli
As noted above, elastic moduli determined for the anterior (CXL-treated) corneal layer can be used to compute effective corneal mechanical moduli, with calculated using a simple rule of mixture, whereas required the inverse rule of mixture. It was interesting for us to check if this model could also describe guided wave behavior in the partially crosslinked cornea when considered as an effective homogeneous material. Using our analytical model, we studied the accuracy of this averaging method to predict effective 'guided wave' moduli, i.e. a pair of and that would fit the analytical -mode dispersion curve of the N-layer case as if the -mode dispersion curve is obtained for a single layer with averaged properties computed from Eqs. (4) and (5). The results are presented in Fig. 6 for a material with a total thickness , made of two layers with and , and and . We studied different repartitions of the layers with either 30% (Fig. 6a), 50% (Fig. 6b) or 70% (Fig. 6c) of the total thickness for the anterior layer respectively. The -mode computed with the effective mechanical moduli using the mixture rules described above does not match the exact analytical solution computed using the individual mechanical moduli of the layers, i.e. using the 2-layered model directly. The difference is especially pronounced in the high-frequency range. Thus, the mixture rules describing effective mechanical moduli in the cornea (or any other layered medium) cannot describe guided wave behavior in the effective medium. This finding is very important and shows that the model of an effective medium cannot be used to describe guided wave behavior. Therefore, reconstruction of effective mechanical moduli from OCE measurements in the cornea should be performed with care. The second question is which mixture model would best describe guided wave behavior in a partially CXL treaded cornea. This is an open and non-trivial question that is outside the scope of this paper. However, as shown in Fig. 6, the use of a simple direct mixture rule for both in-and -out-of-plane moduli (Eq. (5)) produces a dispersion curve closely matching the 2-layer model. In particular, equal thicknesses layers (50/50 split ratio, Fig. 6b) show a near-perfect match. The larger the difference between the layer thicknesses, the larger the difference between the solutions, although the difference remains small. We simulated many different cases, varying thicknesses, moduli, and the number of layers, and can conclude that the mixture rule described by Eq. (5) for both moduli gives a dispersion solution very close to the exact solution. We can treat this observation as an experimental fact but, unfortunately, we do not have its rigorous analytical confirmation. Based on the observation above, we also reconstructed elastic moduli in the anterior (CXL treated) layer of the cornea using Eq. (5) as a mixture rule for both in-and out-of-plane moduli and obtained and , which is in overall good agreement with the directly (using 2-layer model) measured post-CXL effective moduli and . Again, we do not have a mathematically rigorous solution for the 'effective' guided wave behavior but can confirm that both theoretical simulations and experimental results suggest a simple mixture rule for both mechanical moduli to describe the effective guidance of mechanical waves in multi-layered nearly-incompressible media in general, and in the cornea in particular.

IV. DISCUSSION AND CONCLUSIONS
In this study we combined structural OCT with dynamic AT-OCE to assess the penetration depth of CXL treatment in the cornea. Analyzing the brightness of structural OCT images and the rate of image decorrelation between consecutive B-scans, we can conclude that there is a sharp transition between CXL and untreated cornea layers. This finding allowed us to measure the thickness of the CXL treated layer and suggested a model of a two-layer medium to reconstruct both in-and out-of-plane elastic moduli in this layer. Using AμT-OCE, we tracked guided wave propagation in a cornea before and after CXL. Using an appropriate NITI model (in each layer), we quantified the localized stiffening of cornea with CXL. Because the -mode occupies the whole thickness, it carries averaged information about material stiffness (i.e., averaged over its two parts). We proposed two ways to determine the effective moduli in the treated cornea. In the first case, the 2D-spectrum in the post-CXL cornea was fitted with the two-layer model with known thicknesses of each layer and known moduli of the untreated part (obtained from OCE measurements in the untreated cornea). The effective mechanical moduli of the entire cornea can be then calculated using Eqs. (4) and (5) respectively for out-of-and in-plane moduli. The second way assumed the effective guided wave The total thickness is , but the distribution of the layers varies so that: a) the top layer takes 30% of the total thickness and the bottom one 70%; b) both top and bottom layers take 50% of the total thickness; c) the top layer takes 70% of the total thickness and the bottom one 30%. behavior in the cornea, where the 2D-spectrum post-CXL is fitted with a single layer model. We found that the mixture rule for layered materials in continuous mechanics (see Eqs. (4) and (5)) works differently for in-and out-of-plane moduli. From these expressions the effective moduli of partially CXL-treated cornea can be computed, which may be important in assessing CXL surgery and predicting its outcome. However, the same mixture rules cannot be applied to quantify effective guided wave behavior in the layered medium. Both numerical simulations and ex vivo experiments in the human cornea sample suggest a simple mixture rule (Eq. (5)) for both corneal moduli, but we do not have a rigorous analytical confirmation of this observation. One of the limitations of the study is that it assumes the post-CXL cornea consists of two homogeneous layers. This model was found reasonable in several previous studies. In fact, penetration of the riboflavin/dextran solution into the cornea and the unequal UVA irradiation with depth might suggest a more gradual transition of stiffness between layers. A multilayer model following the approach described here can be used to compute the -mode dispersion curve for an arbitrary number of layers, their stiffnesses and thicknesses and, thus, quantify CXL treatment outcomes with better precision. Recent results suggest that reverberant OCE can reconstruct the stiffness depth-dependent variation [29]. It would be interesting to compare our method with reverberant OCE, which will be the object for our future studies. Note, however, that reverberant OCE is not currently feasible in vivo because it uses contact vibrators. This is why guided wave-based OCE is still the only method capable of in vivo non-contact measurements of corneal anisotropic elasticity, ultimately to be spatially resolved. Finally, we have shown that phase-sensitive OCT combined with AT wave excitation can accurately measure both the structure of human corneas and the depth-dependence of moduli due to crosslinking. These findings are essential for building personalized models of corneal deformation following CXL, and thus better adapt crosslinking therapy for clinical use and predict its outcomes. However, further experiments on a larger group of samples are required to generalize the present results.